//
//  matrix.hpp
//  utilities
//
//  Created by Sam Jaffe on 8/15/13.
//  Copyright (c) 2013 Sam Jaffe. All rights reserved.
//

#pragma once

#include "base.hpp"
#include "matrix_helper.hpp"
#include "matrix_R_1.hpp"
#include "matrix_1_C.hpp"
#include "matrix_1_1.hpp"
#include <array>

namespace math {
    namespace linalg {
        
        template <size_t _N_row, size_t _N_col, typename T>
        class matrix : public _matrix<_N_row, _N_col, T, std::array<T, _N_col>> {
        private:
            typedef matrix<_N_row, _N_col, T> _Self;
            typedef _matrix<_N_row, _N_col, T, std::array<T, _N_col>> _Base;
            typedef matrix_column<_N_row, _N_col, const T, const matrix> _c_column;
        public:
            typedef typename _Base::value_type value_type;
            typedef typename _Base::reference reference;
            typedef typename _Base::const_reference const_reference;
            typedef typename _Base::size_type size_type;
            typedef matrix_column<_N_row, _N_col, T, matrix> column_type;
            typedef const matrix_column<_N_row, _N_col, T, matrix> const_column;
            typedef typename _Base::row_type row_type;
            typedef typename _Base::const_row const_row;
            typedef typename _Base::row_vector row_vector;
            typedef typename _Base::column_vector column_vector;
            typedef typename _Base::row_matrix row_matrix;
            typedef typename _Base::column_matrix column_matrix;
        private:
            std::array<std::array<T, _N_col>, _N_row> data;
            
            virtual const _c_column dummy_col(const size_type index) const {
                return _c_column(*this, index);
            }
            
            virtual const_row dummy_row(const size_type index) const {
                return this->operator[](index);
            }
        public:
            constexpr matrix() {
                
            }
            matrix(const _Self& other) = default;
            matrix& operator = (const _Self& other) = default;
            matrix(_Self&& other) = default;
            matrix& operator = (_Self&& other) = default;
            
            matrix(const std::array<std::array<T, _N_col>, _N_row>& other) :
            data(other) {
                
            }
            
            virtual ~matrix() {}
            
            _Self operator + (const_reference c) const {
                return _Base::add(c, _Self{});
            }
            
            _Self operator - (const_reference c) const {
                return _Base::sub(c, matrix{});
            }
            
            _Self operator * (const_reference c) const {
                return _Base::mult(c, matrix{});
            }
            
            _Self operator / (const_reference c) const {
                return _Base::div(c, matrix{});
            }
            
            _Self operator - () const {
                return _Base::neg(matrix{});
            }
            
            /**
             * Ex: D1 = 3, D2 = 2, D3 = 1
             * [ 1  0 ]   [ 1 ]   [ 1 ]
             * [ 0  1 ] * [ 2 ] = [ 2 ]
             * [ 1  1 ]           [ 3 ]
             */
            template <size_t _N_col2>
            matrix<_N_row, _N_col2, T>
            operator * (const matrix<_N_col, _N_col2, T>& u) const {
                matrix<_N_row, _N_col2, T> C;
                for (size_type j = 0; j < _N_col2; j++) {
                    for (size_type i = 0; i < _N_row; i++) {
                        C[j][i] = u.dummy_col(j).dot(this->dummy_row(i));
                    }
                }
                return C;
            }
            
            virtual _Self operator + (const _Self& B) const {
                _Self C;
                for (size_type j = 0; j < _N_col; j++) {
                    for (size_type i = 0; i < _N_row; i++) {
                        C[j][i] = this->operator[](j)[i] + B[j][i];
                    }
                }
                return C;
            }
            
            virtual _Self operator - (const _Self& B) const {
                _Self C;
                for (size_type j = 0; j < _N_col; j++) {
                    for (size_type i = 0; i < _N_row; i++) {
                        C[j][i] = this->operator[](j)[i] - B[j][i];
                    }
                }
                return C;
            }
            
            matrix<_N_col, _N_row, T> transpose() const {
                return _Base::_transpose(matrix<_N_col, _N_row, T>{});
            }
            
            virtual row_type operator[](const size_type index) {
                return data[index];
            }
            
            virtual const_row operator[](const size_type index) const {
                return data[index];
            }
            
            virtual row_type row(const size_type index) {
                return this->operator[](index);
            }
            
            virtual column_type col(const size_type index) {
                return column_type(*this, index);
            }
            
        };
        
    }
}
